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ABSTRACT 

The nonlinear dynamics of outflows driven by magnetic explosion on the surface of a compact star is investi- 
gated through special relativistic magnetohydrodynamic simulations. We adopt, as the initial equilibrium state, 
a spherical stellar object embedded in hydrostatic plasma which has a density p(r) oc r~ a and is threaded by a 
dipole magnetic field. The injection of magnetic energy at the surface of compact star breaks the equilibrium 
and triggers a two-component outflow. At the early evolutionary stage, the magnetic pressure increases rapidly 
around the stellar surface, initiating a magnetically driven outflow. A strong forward shock driven outflow is 
then excited. The expansion velocity of the magnetically driven outflow is characterized by the Alfven veloc- 
ity on the stellar surface, and follows a simple scaling relation v mag oc v^l 1 . When the initial density profile 
declines steeply with radius, the strong shock is accelerated self-similarly to relativistic velocity ahead of the 
magnetically driven component. We find that it evolves according to a self-similar relation T s h oc r s h, where T s h 
is the Lorentz factor of the plasma measured at the shock surface r s h. Purely hydrodynamic process would be 
responsible for the acceleration mechanism of the shock driven outflow. Our two-component outflow model, 
which is the natural outcome of the magnetic explosion, can provide a better understanding of the magnetic 
active phenomena on various magnetized compact stars. 

Subject headings: MHD, relativistic processes, stars: neutron - stars: winds, outflows, methods: numerical 



1. INTRODUCTION 

Plasma outflow from gravitationally bounded stellar objects 
is a universal phenomenon in astrophysics which is known 
to occur on a range of different spatial and energetic scales. 
Gaseous material containing synthesized elements and energy 
are supplied from the stellar source into the surrounding envi- 
ronment by these outflows. It is important to study the driving 
mechanism of the plasma outflow not only for the dynamical 
evolution of the stellar object itself, but also for the chemical 
and dynamical evolution of interstellar medium of galaxies. 

The pioneering theory of stellar wind, which is the proto- 
type of the astrophysical outflow, was developed by Parker 
(1958). It describes steadily accelerated plasma flow from a 
high-temperature corona into interstellar space. Weber & De- 
vis (1967) extended Parker's wind theory to include magnetic 
processes by modeling plasma outflows along open magnetic 
field line. In this model, the accelerated stellar wind collides 
with the interstellar medium and finally forms the shock wave. 

The neutron star is known to have a relativistic plasma out- 
flow, called as pulsar wind. The magnetosphere of rotating 
magnetized neutron star is filled with charged particles (Gol- 
dreich & Julian 1969; Sturrock 1971; Ruderman & Sutherland 
1975). These charged particles flow out steadily along open 
magnetic field lines through the light cylinder where the co- 
rotation velocity coincides with the light speed. The plasma is 
accelerated to relativistic velocity by converting the rotational 
energy of the pulsar itself. The pulsar wind also interacts with 
the surrounding dense gas ejected from the progenitor stars 
and causes a termination shock (Kennel & Coroniti 1984). 
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Non-steady outflows exist universally in astrophysical sys- 
tems. These should be seen as a more general class of out- 
flows than the steady stellar wind. Of these non-steady flows, 
astrophysical jets are the most studied. Astrophysical jets are 
collimated bipolar outflows from central objects. It is be- 
lieved that magneto-centrifugal and magnetic pressure forces 
are the most promising power for driving the collimation flow 
(Blandford & Payne 1982; Uchida & Shibata 1985; Shibata 
& Uchidal986). The magnetized jet from young stellar ob- 
jects have a central role in prompting the collapse of molec- 
ular clouds and thus is essential for star formation (Machida 
et al. 2008). The super-massive black holes in active galactic 
nuclei also powers bipolar jets (Begelman et al. 1984; Bridle 
& Perley 1984). The most energetic event in the universe is 
the gamma-ray burst, and which is believed to be energized by 
ultra-relativistic bipolar flows formed when a massive stellar 
core collapses (Woosley 1993; MacFadyen & Woosley 1999). 

The non-steady mass ejection event associated with solar 
flares is called "coronal mass ejection (CME)" (Low 1996). 
The emergence of a twisted magnetic flux tube from below 
the solar surface, which supplies free energy to the system, 
triggers flares by magnetic reconnection (Ishii et al. 1998; 
Kurokawa et al. 2002; Magara 2003) and drives the associated 
CMEs (e.g., Mikic & Linker 1994; Shibata et al. 1995). 

We can observe non-steady mass ejection phenomena anal- 
ogous to the solar flare/CME activity in magnetar systems. 
The magnetar is a ultra-strongly magnetized neutron star with 
magnetic field of 10 I4 -10 15 G (Duncan & Thompson 1992, 
Thompson & Duncan 1995; Harding & Lai 2006; Mereghetti 
2008). The giant flare observed from the magnetar is surpris- 
ingly energetic ( 1 44 — 1 46 ergs) and is believed to be pow- 
ered by releasing the enormous magnetic energy stored in the 
magnetar itself (Thompson & Duncan 1995; Lyutikov 2006; 
Masada et al. 2010). We note that the magnetar giant flare 
is also associated with an expanding ejecta (Cameron et al. 
2005; Gaensler et al. 2005; Taylor et al. 2005). 

While the energetic origin of the outflow is different in the 
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various astrophysical systems, there is a characteristic com- 
mon to all outflow phenomena. The magnetohydrodynamic 
(MHD) processes play a role in powering and regulating the 
outflow. It has been a central issue in the study of the astro- 
physical outflow to specify the MHD process which controls 
the outflow dynamics (Hayashi et al. 1996; Beskin & Nokh- 
rina 2006; Spruit 2010). Nevertheless, the physical property 
and the acceleration mechanism of the outflow powered by 
the magnetic energy released from the compact stars have not 
yet been studied sufficiently (e.g., Komissarov & Lyubarsky 
2004; Komissarov 2006; Takahashi et al. 2009). 

This paper focuses on the outflow powered by a magnetic 
explosion on a compact star, such as neutron star. We inves- 
tigate its nonlinear dynamics using 2.5-dimensional special 
relativistic MHD simulations. The purpose of this paper is to 
reveal the characteristics of the outflow which expands from 
the surface of the compact star into the surrounding medium. 
The relativistic effect on the physical properties of the outflow 
is a special interest of this work because such outflow should 
be easily accelerated to relativistic velocity. 

This paper opens with descriptions of the numerical model 
in § 2. The numerical results of our relativistic MHD simu- 
lation are presented in § 3. The self-similar property of the 
outflow is studied in more detail using one-dimensional hy- 
drodynamic simulation in § 4. In § 5, the characteristics and 
potential application of our outflow model are discussed. Fi- 
nally, we summarize our findings in § 6. 

2. NUMERICAL MODEL 

2.1. Governing Equations 

We solve the special relativistic magnetohydrodynamic 
(SRMHD) equations in an axisymmetric spherical polar co- 
ordinate system (r, 9, </>)■ The fluid velocity v, electric field 
E and magnetic field B in the system have three components 
while their derivatives in the (^-direction are assumed to be 
zero. We adopt, as the equation of state for the surrounding 
medium, the ideal gas law with a ratio of the specific heats 
7 = 4/3. For simplicity all the dissipative and radiative effects 
are ignored. The basic equations are then 

dD 
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Here T = [1 -(v/c) 2 ] '/ 2 is the Lorentz factor, and e = pc 2 + 
P/( j- 1) is the internal energy. In Newtonian gravity g, we do 



not take account of the relativistic correction on the gravita- 
tional force because it should be negligible in the neutron star 
system. D, R and e denote the density, the momentum den- 
sity and the energy density measured in the laboratory frame 
respectively. The other symbols have their usual meanings. 

We use the HLL scheme (Harten et al. 1983) to solve the 
SRMHD equations in the same manner as Leismann et al. 
(2005). The numerical flux across the cell interface is eval- 
uated using the characteristic velocity of the fast magneto - 
sonic wave obtained from the physical variables in left and 
right adjacent cells of the interface. The primitive variables 
are calculated from the conservative variables following the 
method of Del Zanna et al. (2003). We use a MUSCL-type 
interpolation method to attain second order accuracy in space, 
and Constrained Transport method to guarantee a divergence 
free magnetic field (Evans & Hawley 1988). 

2.2. Initial Setting of Numerical Model 

For our initial conditions, we consider a hydrostatic gas sur- 
rounding a central compact star. A dipole magnetic field an- 
chored in the central star threads the external gaseous plasma, 
and is given by 



B r = 



2B Rjcos6 



Be- 



BqR* sinO 



B 6 =0. 



(9) 

(10) 
(11) 



where Bq and are the field strength at the equatorial plane 
on the stellar surface and the radius of the central star respec- 
tively. We assume that the temperature of the surrounding 
medium is inversely proportional to the distance from the cen- 
troid of the star 



(12) 



r = r (- 



where Tq is the reference temperature of the surrounding gas 
on the surface of the central star. Then the spatial distribution 
of the gas density and pressure in hydrostatic equilibrium are 



P = Po 



P = Pn 



R* 



where 
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PoR* 
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(13) 
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Here po and Pq are the reference density and pressure of the 
surrounding gas on the central star, and M* is the central stel- 
lar mass. The parameter a denotes the power-law index of the 
initial density profile. 
The convectively stable condition for the surrounding gas is 



ds 

dr >0 



(16) 



where s is the entropy (Kippenhahn & Weigert 1990). When 
adopting equations (12)— (15), this can be reduced to 



a > 



1 



7-1 



= 3 . 



(17) 
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We should set the parameter a to be larger than 3 in our cal- 
culations for maintaining the surrounding plasma being con- 
vectively stable. On the other hand, models with a > 7.25 can 
not be solved using our simulation code because the physical 
variables become too small to be resolved in the calculation 
domain. The parameter a surveyed in the following is thus 
restricted to the range between 3.0 and 7.25. 

The ratio between the mass and radius of the central star 
is a fundamental parameter which characterizes the system, 
and is fixed as M*//?* = 2.2 x 10 27 g/cm. This is the typical 
value for the neutron star. All the numerical models we exam- 
ined here are controlled by two non-dimensional parameters, 
which are the power-law index a of the hydrostatic surround- 
ing medium and the plasma beta on the equatorial surface of 
the central star /3 = Po/(BI/8tt). 

As we show in § 3.2, the temporal evolution of the magneti- 
cally driven outflow is controlled by the initial Alfven velocity 
(va). The larger value of initial plasma beta (oc v A 2 ) provides 
the slower evolution of the outflow, which is not appropriate 
for extensive parameter survey. In addition, the conservative 
MHD scheme used in our code is not good at solving the prob- 
lem with small values of the initial plasma beta. Hence we 
survey the /^-dependence of the outflow properties within the 
range 3.2 < A) < 3.2 x 10 3 . 

The relativistic representations of the sound speed and 
Alfven velocity v s s and v* are, at the equatorial surface of 
the central star, 
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Here v s ,o and va.o are the sound speed and Alfven velocity at 
the equatorial surface of the central star in the non-relativistic 
limit. Since the parameter a varies from 3.0 to 7.25 and (3q 
varies from 3.2 to 3.2 x 10 3 in our models, the relativistic rep- 
resentations of the sound and Alfven velocities can be well 
approximated by their non-relativistic forms as stated above. 
Using description (fTBI l, the sound and Alfven velocities are 



v s0 / c = 0.48(a+iy 



-1/2 



rWr = 0.32<n + ir' 2 ( || 



-1/2 



(21) 
(22) 



In the following, we mainly use the non-relativistic form of 
the sound and Alfven velocity for convenience. The normal- 
ization units in length, velocity and time are chosen for all the 
models as radius of the central star /?*, light speed c, and the 
light crossing time t=R*/c. 

2.3. Boundary Condition and Numerical Grid 

The stress-free condition is applied to the outer radial 
boundary located at r = 200/?*. At the inner radial boundary 
(r = /?*), we impose the condition that the physical variables 
except and B$ are fixed to be time-independent, that is, 



d/dt\ r= R t = 0. The radial and latitudinal velocities of the gas 
v r and vg vanish at the stellar surface. The gas density, the 
gas pressure and the poloidal field strength at the stellar sur- 
face are consistently determined from equations (9)-(14). The 
stress-free condition is adopted at the inner boundary only for 
the toroidal magnetic field B^. 

The shearing motion at the inner boundary regulates the en- 
ergy injection into the calculation domain. On the equatorial 
surface, we impose a shear flow with the longitudinal velocity 



v*(**, 6, t) = Av (f )6 exp[( 1 - 6 4 )/4] 



(23) 



where 9 = (8 - 7r/2)/A6» peak (see Mikic & Linker 1994). 
A6> pea k controls the latitudinal width within which the shear- 
ing motion is operated, and is fixed as tt/9 in the following. 

We investigate two different energy injection patterns. In 
both cases, the shearing velocity Av<f,(t) increases linearly 
from zero to 0.1c between t = and lOr and is held at con- 
stant velocity 0.1c until t = 20r. After t = 20r, the velocity of 
the shearing motion is kept constant at 0.1c in case A, and is 
terminated suddenly in case B. 

The amplitude of the shearing velocity is an important pa- 
rameter which controls the outflow dynamics, and should be 
widely surveyed. However, we focus only on the typical case 
with the fixed shearing velocity of 0.1c in this paper. In this 
case, the shearing motion supplies the magnetic energy which 
is roughly comparable to the thermal energy of the system. 
The dependence of the outflow dynamics on the initial field 
strength and the density profile of the surrounding gas is our 
interest in this paper. We will study the impact of the shearing 
velocity on the outflow dynamics in the separate paper. 

The polar angle of the calculation domain is uniformly di- 
vided into A6> = 7r/90. A nonuniform grid is adopted in the 
radial direction. The radial grid size Ar gradually increases 
from 0.0361?* (at r = /?,) to 0.2/?* (at r = 5/?*), and is fixed 
as 0.2/?* outside r = 5/?*. The largest grid size Ar = 0.2/?* is 
empirically determined to correctly capture the the sharp rela- 
tivistic shock which propagates through the outer calculation 
domain. The number of the grid points is set as 1024 in the 
radial and 90 in the direction for all the models. 

2.4. Dynamical Effects of the Boundary Condition 

The shearing motion imposed on the equatorial surface 
of the compact star generates the toroidal magnetic field by 
stretching the initial dipole field. The kinetic energy of the 
shearing motion is then immediately converted to magnetic 
energy. The magnetic energy amplified above the stellar sur- 
face initiates the dynamical evolution of the system. While the 
shearing motion is added artificially to break the initial equi- 
librium of the system, it is physically modeling the supply of 
the magnetic helicity and energy from the stellar interior into 
the magnetosphere, like as the flux emergence in the sun. 

In the solar flare/CME process, the emergence of the tightly 
twisted magnetic flux plays a triggering role in breaking the 
equilibrium of the corona (Ishii et al. 1998; Kurokawa et al. 
2002; Magara 2003). We use the shearing motion instead of 
the flux emergence for the sun as the process modeling energy 
and helicity injection in our calculations, because actual en- 
ergy and helicity injection processes are still veiled in mystery 
in the context of compact stars. 

Since the inner boundary condition in our models is math- 
ematical over-determined (see Bogovalov 1997), there exists 
a discontinuity in the radial velocity at around the stellar sur- 
face in our calculations. Figure QJ depicts the radial profile 
of the radial velocity of the model where a = 3.5 and flo = 3.2 
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FIG. 1 . — Panel (a): the radial profile of the radial velocity for the model where a = 3.5, /3q = 3.2 at the steady state on around the inner boundary region. The 
filled square, triangle and circle represent the profiles at the different polar angles 8 = tt/6, it/3, and tt/2 (along the equatorial plane), respectively. Panel (b): the 
density distribution (log scale) and poloidal velocity vector near the inner boundary on the meridional plane at the steady state. The density is normalized by its 
inner boundary value. The orientation and length of each arrow denote the direction and magnitude of the poloidal velocity vector at each meridian point. 




tlx 

FIG. 2. — The temporal evolution of the total mass contained in the cal- 
culation domain for the model where a = 3.5, /3q = 3.2. The vertical axis is 
normalized by the initial total mass of the plasma distributed in the entire 
domain. 

as a reference. Here we focus on the inner boundary region 
(1 < r/R* <1.2) and the time-steady state of t = 350t. The 
filled square, triangle and circle represent the radial velocity 
profiles at different polar angles 6 = tt/6, ir/3, and ir/2 re- 
spectively. This addresses that there exists a velocity discon- 
tinuity near the inner boundary. In addition, the discontinuity 
at around the high latitude stellar surface has a larger ampli- 
tude than that around the equatorial region where the outflow 
is mainly accelerated in our models as shown in § 3. 

Figure QJ shows the meridional distribution of the density 
and poloidal velocity near the inner boundary at the steady 
state of the model same as Figure QJ. The orientation and 
length of each arrow denote the direction and magnitude of 
the poloidal velocity vector at each meridian point. The grey 
contour is a logarithmic representation of the density normal- 
ized by its inner boundary value. The radial flow across the 
inner boundary is induced especially at around the high lat- 
itude stellar surface although the radial velocity is set to be 
zero at the inner boundary. This induced outgoing radial flow 
supplies the plasma into the calculation domain. 

The temporal evolution of the total mass contained in the 



calculation domain is presented in Figure|2]in the same model. 
The vertical axis is normalized by the initial total mass of the 
plasma distributed in the entire domain. The total mass in 
the system increases with time due to the mass inflow from 
the inner boundary region. The mass injection rate is roughly 
10~ 2 Mini/r at the initial evolutionary stage until IOOt, but 
10" 3 Mi n i/T at the steady state after 200t. 

Since the density of the plasma near the stellar surface is 
maintained at close to the initial value by the mass inflow 
across the inner boundary as shown in Figure QJ, the Alfven 
and sound velocity (oc p~ 1//2 ) are not relativistic there and an 
order of 0.1c at the steady state. The radial velocity is less 
than 0.01c when focusing on the equatorial region where the 
magnetic energy is strongly amplified and the outflow is ener- 
gized, while its latitudinally-averaged value is roughly 0.05c. 

3. RESULTS 
3.1. Temporal evolution of a fiducial model 

We show the dynamical evolution of the outflow in a fidu- 
cial model (a = 3.5 and (3q = 3.2). The case with continuous 
shearing motion (case A) is investigated here. 

Figure [3] demonstrates the temporal evolution of the ex- 
panding outflow. The density and the magnetic field lines 
projected on the meridional plane are shown at the time (a) 
t = 0, (b) 350r and (c) 700t respectively. Note that the grey 
contour is a logarithmic representation of the density normal- 
ized by its initial value. 

The magnetic pressure resulted from the toroidal magnetic 
field, which is generated by the shearing motion around the 
stellar equatorial surface, initiates and accelerates the quasi- 
spherical outflow. The driven outflow is strongly magnetized, 
in which the magnetic pressure composed mainly of the con- 
tribution for the generated toroidal field dominates the gas 
pressure, and carries dense plasma away from the central star. 

In the equatorial region where the magnetized outflow is en- 
ergized, it is natural from the physical perspective to consider 
that the magnetization parameter (a = B 2 / 'AttT 2 pc 2 ) grows 
due to both the mass ejection and the continuous injection 
of the magnetic energy. However, the magnetization parame- 
ter near the stellar surface is maintained at a lower value than 
unity during the computing time as shown in Figure |4] 
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FIG. 3 . — The time evolution of the density distribution (log scale) and magnetic field lines in the meridional plane for the model where a = 3.5 and 0q = 3.2. 
The density is normalized by its initial value: p(r, 9, t)j p(r, 9,0). The left, middle and right panels are corresponding to those in t/r = 0, 350, and 700 respectively. 
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panel. 



6 



Matsumoto et al. 



10 



10 



Q_ 10" 



Q. 



10 



10" 



10 



-10 



1 1 1 1 1 1 1 1 1 1 1 w- 

i (n\ 
1 w 


- 


- (=350r 


- 

r = 700r' 












50 



100 

r/R. 



150 200 



10 



oc 



~o 10 



oc 



K 



10 
10" 
10" 
10 

0.4 

0.3 

0,2 

0.1 

0.0 



■10 



(b) 



t = 700x' 



r = 350x 



/=0 



50 



100 

r!R, 



150 200 





(c) 







~f\ 

1 V ^ _ / 


/ = 700t 




J *= 350x 






1 


t = 





50 



100 



150 



200 



FIG. 5. — Panels (a) the density, (b) the total pressure (thermal + magnetic 
pressures), and (c) the radial velocity on equatorial plane as a function of the 
normalized radius r/fi* where a = 3.5 and 0g = 3.2. Dotted, dashed, and 
solid lines are corresponding to those in t/r =0, 350 and 700, respectively. 



In this figure, the spatial distribution of the magnetization 
parameter (left) and density (right) are demonstrated by the 
logarithmic grey scale on the meridian plane. Different panels 
are corresponding to the different evolutionary phases (a) t = 
0, (b) 100t, (c) 300t and (d) 700t respectively. The region 
occupied by the central star is filled by black and white in 
the left and right parts of each panel. The magnetization a 
near the stellar surface is less than unity. The mass inflow 



from the inner boundary region plays a role in reducing the 
magnetization parameter to less than unity. 

Although the magnetic energy is strongly amplified and the 
outflow carries a large portion of the gas around the equatorial 
region, the rapid amplification of the magnetization parameter 
is prevented because the dense plasma maintained above the 
high-latitude stellar surface is continuously supplied into the 
equatorial region. This is the reason why the magnetization 
parameter a in the near surface region is less than unity at all 
the latitude during the computing time. 

Figure [5] denotes the radial distributions of the (a) density, 
(b) total pressure (thermal + magnetic pressures) and (c) ra- 
dial velocity along the equatorial plane (9 = tt/2). The dot- 
ted, dashed and solid lines represent the cases t = 0, 350r and 
700t of Figure [3] respectively. Note that the vertical axis is 
normalized by its reference value in each plot. 

It is found that a strong forward shock is formed and 
expands supersonically into the surrounding medium. The 
shock powered by the magnetized outflow is accelerated to 
sub-relativistic velocity ~ 0.3c. Behind the forward shock 
surface (around 50/?* when t = 700r), a reverse shock is ex- 
cited and propagates through the shocked medium. 

The time evolution of the gas and magnetic pressures are 
shown in Figure [6] The left column is of the case being dis- 
cussed in this section. The top, middle and bottom panels 
correspond to phases (a) t = 0, (b) 350r, and (c) 700r in the 
fiducial model. In each panel, the solid line indicates the gas 
pressure, dotted and dashed lines are the magnetic pressure 
due to the poloidal and toroidal magnetic field respectively. 
The normalization of each curve is the initial magnetic pres- 
sure at the equatorial surface of the central star, that is Bq 2 /8it. 

In Figure [6}:, we find a tangential discontinuity created be- 
hind the propagating forward shock. The magnetic pressure 
provided by the toroidal magnetic field becomes predominant 
behind the discontinuity. In contrast, the gas pressure domi- 
nates the magnetic pressure between the discontinuity and the 
shock surface throughout the evolution of the system, except 
for a transition region near the discontinuity. 

3.2. Magnetically Driven Outflow 

The physical properties of the magnetically driven outflow 
are examined through changing the a-parameter that controls 
the density profile. We anticipate that the magnetically driven 
outflow velocity v mag is characterized by the strength of the 
initial dipole field which is the seed of the predominantly 
toroidal field behind the tangential discontinuity. The rela- 
tion between the velocity of the magnetically driven outflow 
and the initial field strength is elucidated here. 

Figure [7] shows the outflow velocity v mag in relation to the 
initial Alfven velocity at the equatorial surface of the central 
star (va.o) for the models a = 3.5, 4.5, 5.5, 6.5, and 7.0 respec- 
tively. The initial Alfven velocity is the representative of the 
seed field strength and varies from 1 .1 x 10~ 3 c to 0.15c when 
changing the size of a and (3. The velocity is measured when 
the head of the magnetically driven outflow reaches r = 70/?*. 

The velocity of the magnetically driven outflow increases 
with the initial Alfven velocity. The logarithmic fitting of 
the data provides the power law dependence between them 
as summarized in Table 1, where a = dlnv mag /dlnvA.o de- 
notes the power law index. The index a is nearly constant for 
all the models, that is a ~ 0.5. The dynamical evolution of 
the magnetically driven outflow thus follows a simple scaling 
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FIG. 6. — Left column: temporal evolution of the radial distribution of the gas and magnetic pressures where a = 3.5 and /3o = 3.2 at the times (a) t = 0, (b) 
t = 350t and (c) t = 700r respectively. Right column: the model a = 7.25 and /3q = 3.2. Note that all pressure components are measured on the equatorial plane. 
The magnetic pressure is divided into two parts in these panels, one is the contribution from toroidal magnetic fields B,, and the another is that from the poloidal 



8 



Matsumoto et al. 



10' 



10" ■ 



10" 





"i 


— ■ — ■ — ■ — ■ - - - 1 ■ — - — - — 














a 










' * 




K 








A/ 


* 


a 3 .5 


♦ 






cc = 4.5 


• 


♦ 




a = 5.5 








a = 6.5 


□ 






a = 7.0 


X 






Vmag « VA.0 





10 



10" 



10" 



10' 



VVui/t' 



Fig. 7. — The relation between the Alfven velocity at the equatorial surface 
and the velocity of the tangential discontinuity measured at r = 70i?„ when 
a = 3.5, 4.5, 5.5, 6.5, and 7.0. There is a simple scaling relation between 
these two values, v m ag <x va CI - 



TABLE 1 

Indices of a simple relation between the magnetically 
driven outflow velocity at tangential discontinuity surface 
and the initial alfven velocity at the equatorial surface of 
the central star when the power law index a of the initial 
density profile varies. the relation is v mag oc v a o ct . 



a 3.5 


4.5 


5.5 


6.5 


7.0 


cr 0.51 


0.54 


0.45 


0.44 


0.48 



To draw a qualitative physical picture which is accountable 
for this scaling relation, we focus on two balancing equations. 
For simplicity, we reduce the governing equations (l)-(8) to 
their non-relativistic forms because the outflow velocity of in- 
terest here is sufficiently slower than the light speed. The nu- 
merical results address that in the steady state, the strength of 
the toroidal field behind the tangential discontinuity does not 
change significantly in time. Hence the radial advective loss 
of the toroidal field should approximately counterbalance the 
generation of the toroidal field by shearing motion at around 
the stellar surface. The azimuthal component of the induction 
equation can be simplified, at the steady state, 



dt 



■ = [V x (y x 

B()Av'0 ^ fi<ft,aropVr,adv 



(24) 



where B$ Amp is the typical strength of the toroidal field in the 
region where it is strongly amplified and v,.. a dv is the mean ra- 
dial advection velocity there. The first term on the right hand 
side of equation d24T i represents the generation of the toroidal 
field by the shearing motion within typical latitudinal width L, 
which is replaced using L ~ A6> pea k in the following. The 
second term shows the advective loss of the toroidal field with 
the mean radial velocity v,. a( jv This yields the mean toroidal 
field in the steady state, 



B 



<5>,amp - 



l 



Aft 



peak 



^r.adv / 



5„ 



(25) 



Our numerical results suggest that the magnetically driven 
outflow is mainly accelerated in the near-surface region where 



the toroidal field is strongly amplified as shown in Figure |6j; 
and [7J;. The velocity of the magnetized outflow does not 
change drastically in the region beyond the radius r ~ 10/?*. 
This would address that large portion of the kinetic energy of 
the outflow is gained through the conversion of the magnetic 
energy in the near-surface region of the compact star. 

In the steady state, the total energy of the magnetically 
driven outflow should be conserved along the streamline, 
which is almost radial in our models. The specific kinetic en- 
ergy of the outflow at an arbitrary radius, at which the energy 
conversion is almost terminated, is thus given by integrating 
the magnetic pressure gradient force along the radial stream- 
line from the stellar surface to the arbitrary radius r, 

2' 



1 



M r )=- 



i 



(26) 



(see Appendix for details). Since the magnetized outflow is 
not accelerated dramatically after the termination of the en- 
ergy conversion, its mean radial velocity Vmag can be substi- 
tuted for the local value v r (r) in equation d26b . This energy 
balance equation yields the typical velocity of the magnetized 
outflow as a function of the mean toroidal field B^.amp and 
the mean density p amp in the region where the toroidal field is 
strongly amplified and the outflow gains the kinetic energy, 

B ' (27) 



'0,amp 



\/47rp am p 

As the outflow velocity does not change drastically in the re- 
gion beyond the near-surface of the central star, the mean 
advection velocity v r ,adv appeared in equation $25[ should 
be comparable to the typical velocity of the magnetically 
driven outflow, that is v r . a dv — 
tions d27l > with d25l >. we can provide 

^Av// 2 Aft 



Then, by combining equa- 



va.o 



1/2 
peak 



Pamp 
P0 



-1/4 



(28) 



where the parameters va,o> Av^,, and A6> pea k are initially given 
and remained to be constants. The relation (f28) obtained us- 
ing an order of magnitude analysis here is shown by the solid 
line in Figure [7] The numerical results are well reproduced 
by this analytic scaling. It is noted again that two balancing 
equations, which are reduced from the induction and energy 
equations, account for the scaling relation. 

The physical property of the magnetized outflow studied 
here is different from that found in the similar work on so- 
lar CME by Mikic & Linker (1994). This would be because, 
in their work, the shearing motion which models the energy 
injection in the solar flare/CME process is much slower than 
the Alfven velocity in comparison with our models. The mag- 
netic energy is gradually stored above the solar surface and is 
released abruptly in their CME models. In contrast, in our 
models, the magnetic energy is immediately converted to the 
kinetic energy without being stored. Although there is less 
information about the injection mechanism and supply rate 
of the magnetic energy in the compact star, we expect that 
the rapid injection of the magnetic energy into the magneto- 
sphere might be the origin of giant flares and the associated 
mass ejection in the magnetar system. We apply our numeri- 
cal model to the magnetar system in § 5. 

3.3. Relativistic Self-Similar Shock 

The second component of the outflow is the shock driven 
outflow which proceeds the magnetically driven component. 




FIG. 8. — The relation between the initial density profile of the surround- 
ing medium (p oc r~ Q ) and the plasma velocity at the forward shock surface 
on the equatorial plane. Note that the velocity is measured when the shock 
surface reaches the equatorial radius r = 150R». The vertical axis represents 
Tuo — 1 = [1 — (vi5o/c) 2 ]~'' 2 — 1 where Tiso and V150 are the Lorentz factor 
and the plasma velocity at the shock surface when it reaches the fixed radius 
r=150ft„. 



FIG. 9. — The time evolution of the Lorentz factor of the expanding gas in 
the model with a = 7.25 and /3q = 3.2. The dashed-two dotted, dot-dashed, 
dashed, and dotted curves indicate the models when r/r = 80, 120, 160, and 
200 respectively. The solid line represents the time trajectory of the Lorentz 
factor of the fluid velocity at the shock surface measured in the laboratory 
frame. The scaling relation between the Lorentz factor T s h and spherical 
radius r s ^ of the shock surface in the laboratory frame is oc r s ^. The 
relativistic shock evolves self-similarly. 



We examine how the initial density profile of the surrounding 
gas affects the shock driven outflow. The plasma beta at the 
equatorial surface is fixed as (3q = 3.2 here for simplicity. 

Figure [8] shows the plasma velocity of the shock surface 
measured at the reference point r = 150/?* and 9 = tt/2 as a 
function of the index a which controls the density profile. 
The vertical axis denotes Ti5o- 1, where Fiso is the Lorentz 
factor defined by [1 - (vi5o/c) 2 ]] -1 ' 2 , and V150 is the plasma 
velocity of the shock at the reference point. The value T — 1 
measured in this figure is useful for treating both relativistic 
and non-relativistic flows simultaneously. When the flow ve- 
locity reaches the ultra-relativistic regime, it roughly provides 
the Lorentz factor itself, that is T— 1 ~ T. In contrast, it is the 
indicator of the flow velocity also in the non-relativistic case 
v <C c, that is T— 1 ~ (v/c) /2. It is found that the shocked 
plasma is eventually accelerated to mildly relativistic veloc- 
ity in the model a = 7.25. There is a correlation between the 
shock velocity and the index a. The a-dependence of the 
Lorentz factor changes at a ~ 5. These indicate that the den- 
sity profile of the surrounding gas is key to reveal the acceler- 
ation mechanism of the forward shock in our model. 

The right column in Figure [6] depicts the radial profile of 
pressure components along the equatorial plane at the time 
t = 0, 60t and 180t for the model with the steep density pro- 
file (a = 7.25). The curves in panels (e), (f) and (g) have the 
same meaning as those for the model a = 3.5. The toroidal 
magnetic field generated by the shearing motion initially pow- 
ers the magnetized component of the outflow. The forward 
shock is then excited and propagates through the surrounding 
medium in the same way as the model with a = 3.5. 

There is a clear difference between the models with a = 3.5 
and a = 7.25. That is the radial structure of the poloidal mag- 
netic field. In the model with a = 3.5, the poloidal field is 
accumulated in front of the tangential discontinuity. The ra- 
dial pressure gradient of the accumulated poloidal field seems 
to push the shock outward. This suggests that the shock wave 
is mainly accelerated by the magnetic pressure gradient force 
when the density profile moderately declines {a < 5). 

In contrast, there is no steep gradient in the poloidal mag- 
netic field behind the shock in the model a = 7.25. Rather 



TABLE 2 

Indices of a relation between the Lorentz factor r sh and the 

EQUATORIAL RADIUS r sh AT THE SHOCK SURFACE IN THE LABORATORY 
FRAME WHEN THE POWER LAW INDEX a OF THE INITIAL DENSITY 
PROFILE VARIES. THE RELATION IS F sh - 1 OC r s h S ■ 



a 5.75 


6.00 


6.25 


6.50 


6.75 7.00 


7.25 


8 0.99 


1.03 


1.05 


1.08 


1.10 1.11 


1.10 



the propagation surface of the poloidal field coincides with 
the shock surface. This suggests that the excited shock wave 
propagates with the accumulated poloidal field when the den- 
sity profile declines steeply (a > 5). 

Temporal evolution of the radial velocity profile in the 
model a = 7.25 is presented along the equatorial plane in Fig- 
ure [9] The vertical axis shows the Lorentz factor of the ex- 
panding gas. The horizontal axis is the equatorial radius. The 
dotted, dashed, dash-dotted, and dashed-two dotted curves 
correspond to the cases f/r = 80, 120, 160 and 200 respec- 
tively. The solid curve denotes the time trajectory of the 
Lorentz factor of the plasma velocity measured at the shock 
surface in the laboratory frame. The evolution of the relativis- 
tic shock has self-similar properties. The fitting of numerical 
data provides a scaling relation between the Lorentz factor T s h 
and the equatorial radius r s h at the shock surface 

T s h oc r sh . (29) 

The Lorentz factor of the plasma velocity at the shock surface 
increases linearly with the equatorial radius. 

The other models with steep density gradient of a > 5.5 
commonly yield a relation close to equation d29| >. Table 2 
summarizes the power-law index S when fitting the numerical 
models by a power-law relation r s h- 1 cx r S h S . It is important 
to stress that all the models which follow the scaling relation 
d29l have relativistic velocity, that is, T\so > 2, at the refer- 
ence point. The models with a < 5.5 do not reach relativistic 
velocity at the reference point and not follow the relation (|29l ). 
In § 5, the self-similar solution numerically discovered in our 
work is compared to the other analytic solutions, including 
the Blandford-McKee solution (Blandford & McKee 1976). 
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At the initial stage, the magnetic pressure, gas pressure, and 
gas density decrease with radius according to r" 6 , r~ a ~ l , and 
r~ a respectively [see equations d9i>-([T4l>l. The plasma beta 
of the surrounding gas finally becomes smaller than unity for 
the model with a > 5 at a great distance even if it is larger 
than unity at the stellar surface. Moreover, the Alfven veloc- 
ity of the initial surrounding gas increases with radius only 
for the model with a > 6. These imply that a critical point, 
at which the MHD effect on the outflow dynamics will be al- 
tered, might exist at around a = 5-6. Though the MHD effect 
might lead to some change in dynamics of the outflow, we will 
show that pure hydrodynamic mechanism can be responsible 
for the self-similar property of the shock by using a simplified 
one-dimensional hydrodynamic simulationl in § 4. 

Our numerical models indicate that the two-component out- 
flow would be a natural outcome of the magnetic explosion 
on a compact star. In the early stages, the magnetically driven 
outflow forms as a consequence of the rapidly rising mag- 
netic pressure around the stellar surface. The magnetized 
outflow loaded with a dense circumstellar plasma expands, 
and finally excites a relativistic forward shock which has self- 
similar properties. Our two-component outflow model is ap- 
plied to the magnetar system in § 5.2. 

3.4. Impact of Energy Injection Pattern 

We clarify the impact of the energy injection pattern at the 
inner radial boundary on the evolution dynamics of the shock. 
As stated in § 2, we impose two different kinds of the energy 
injection pattern. These are the continuous injection case with 
the fixed shearing motion (case A) and the impulsive injection 
case in which we stop the shearing motion after t = 20r (case 
B). The plasma beta at the equator is fixed again as f3o = 3.2. 

Figure[lO]depicts the radial distribution of the density, total 
pressure and Lorentz factor of the outflow measured at the 
equatorial plane for both cases when t = 200r. The dashed 
and solid curves represent the cases A and B respectively. We 
choose the initial density profile with a = 7.25 for this figure. 
While there are differences in the structure of the magnetically 
driven outflow far behind the shock surface, we can not find 
any difference among two cases around the forward shock. 

There is little difference in the structure around the shock 
between two cases when a > 5.25. However, for the model 
with a < 5.25, the different injection patterns provide differ- 
ent spatial structures of the shock as in the structures of the 
magnetically driven component. This can be verified by Fig- 
ure [T2] which shows the radial distribution of the density, to- 
tal pressure and velocity in the fiducial model (a = 3.5) when 
t = 700r for the case A and case B. The curves have the same 
meanings as those in Figure [TT] 

We define a time difference value of the radial velocity on 
the equatorial plane Avv as 



10 



Av r (t ,r) = v r (t,r)- v r (t -At,r), 



(30) 



where At denotes the numerical timestep. The time-distance 
diagram for the time difference value Av r (t , r) is presented for 
the models a = 7.25 (Fig.[L2l and a = 3.5 (Fig.QjJi. The ver- 
tical and horizontal axes denote the normalized time and ra- 
dius. For specifying the temporal evolution of the outflow, the 
radial profiles of Av r are plotted at every reference timestep 
on the time-distance plane. The reference timestep is set by 
the light crossing time T = R 1e /c. Panel (a) corresponds to the 
case A and panel (b) is the case B in each figure. 

In the case B, a wave is excited at t = 20r which corre- 
sponds to the termination time of the energy injection while 
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FIG. 10. — The radial structures of the density, total pressure, and Lorentz 
factor of the plasma as a function of the normalized radius r/R t when a = 
7.25, ft) = 3.2, and t = 200t. Case A, which is expressed by the solid line, 
corresponds to the case with continuous energy injection due to the shearing 
motion. Case B, which is depicted by the dashed line, indicates the case with 
impulsive energy injection. The shearing motion is terminated when t = 20t 
in Case B. 



it is not seen in the case A. This wave would thus transmit 
the information related to the energy injection pattern to the 
surrounding gas. We call this wave as "information wave" 
hereafter. The information wave propagates through the sur- 
rounding medium, but never reaches the forward shock in the 
model a = 7.25 (see Figure [12b). This is because the propa- 
gation velocity of the forward shock is roughly light speed c 
while that of the information wave is about 2c/3. 
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(b) a = 3.5 case A: continuous - 
t = 700x case B: impulsive ■ 



50 



100 150 

r!R, 



200 



0.4 



0.3 



0.2 



0.1 



0.0 



1 1 1 T 1 1 1 1 1 1 1 T T 1 1 1 1 1— 

(c) a =3.5 case A: continuous 




I = 700x case B: impulsive 










rJ «.-" rn 






f 
f 







50 



100 150 

riR* 



200 



FIG. 1 1 . — Almost same as Fig 10, but where a = 3.5 and t = 700r. 



Figure [12c is almost same as Figure fTZb . but follows the 
long-term evolution of the time difference value. The excited 
information wave never interacts with the forward shock even 
when we follow its long-term evolution. This indicates that 
the inner boundary condition does not affect the evolution dy- 
namics of the shock once it has been excited when the initial 
density profile is sufficiently steep (a > 5.25). 

It would be natural to consider the effect of the termina- 
tion time of the shearing motion on the shock dynamics in the 
impulsive injection case. There is actually a critical time f cr j t 
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FIG. 12. — Time-distance diagram for the time difference of the ra- 
dial velocity Av r (t,r) on the equatorial plane in the model a = 7.25 where 
Av,(l,r) = v,-(t,r) — v, (t — At,r). Panel (a) corresponds to the continuous in- 
jection case (Case A), and (b) is the impulsive injection case (Case B). A wave 
which transmits the information for the termination of the energy injection is 
excited at t = 20t in the case B. Panel (c) just describes the long-time evolu- 
tion of Av r in case B. The information wave never reaches the propagating 
forward shock. 



beyond which the termination does not impact on the shock 
dynamics, that is t a ^ ~ 10t. When stopping the shearing mo- 
tion at the time before f CI i t , the information wave can catch up 
with the shock and change its dynamics even if the initial den- 
sity profile is sufficiently steep. This is because the shearing 
velocity at the inner boundary Av^(f ) increases linearly from 
zero to 0.1c during < t < 10t (see § 2.3). The energy injec- 
tion rate which determines the initial shock velocity becomes 
small if the shearing motion is terminated before lOr. This 
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FIG. 13. — Almost same as Fig 12, but when a = 3.5. In this model, the 
information wave can reach the propagating forward shock. 



is the reason why there exist the critical time. The shock dy- 
namics does not change even if we stop the shearing motion 
at the time shorter than 20r as long as it continues longer than 
the critical time. 

4. ONE-DIMENSIONAL HYDRODYNAMIC SIMULATION 

We examine the self-similar properties of the shock driven 
relativistic outflow in more detail using one-dimensional (ID) 
special relativistic hydrodynamic (HD) simulations, which 
enables us to study the numerical model in a wide parameter 
range. We specify here that a purely hydrodynamic process 
is responsible for the self-similar properties of the relativistic 
shock discovered in the two-dimensional (2D) MHD models. 

4.1. Numerical Setting 

We solve the ID spherically symmetric system with 1024 
grid points radially using the same grid spacings as those in 
the 2D-MHD model. The MHD effects are fully ignored by 
dropping all the related terms in governing equations (l)-(8). 

The numerical settings adopted in 1D-HD model are the 
same as those adopted in the 2D-MHD model except for the 
magnetic field. We consider a central compact star surround 
by a hydrostatic gas with the temperature, density, and pres- 
sure distributions described by equations (12)— (14). Since the 
sound velocity follows equation d2Tb with the index a vary- 
ing from 3.0 to 7.25, it becomes non-relativistic (~ 0.2c) at 
the stellar surface. The symbols used in this section have the 
same meanings as those in the 2D-MHD model (see § 2.2). 

The stress-free condition is applied to the outer boundary 
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FIG. 14. — The Lorentz factor of the shocked plasma velocity at r = 150R* 
as a function of the density power-law index a. The cross denotes the model 
calculated by one dimensional (ID) hydrodynamic simulations. The filled 
circle indicates the 2D special relativistic MHD model shown in Fig. 8 again 
for the comparison. The meanings of vertical and horizontal axes are same 
as Fig. 8. Note that F lc f is a reference value of the inflow energy flux which 
corresponds to the case with the inflow velocity 0.1c. 



located at r = 2007?,. At the inner boundary (r = R*), we em- 
ploy the condition that the physical variables except velocity 
are fixed. A plasma inflow is imposed on the inner bound- 
ary in 1D-HD models as a substitute for the shearing motion 
which is the energetic origin of the outflow in 2D-MHD mod- 
els. This plasma flow excites the shock-driven outflow in the 
ID models. The relation between the inflow velocity Vi„ and 
the injected energy flux F; n is given by fj n = p§v\j2. We study 
the properties of the shock driven outflow in the models with 
the inflow velocities v m /c= 1/10, 1/30, and 1/100. Their 
dependence on the index a is also our interest in this section. 

4.2. Properties of the Shock Driven Outflow 

The continuous inflow powers a shock driven outflow. The 
cross symbol in Figure [14] shows the plasma velocity of the 
shock surface measured at r = 150/?* as a function of the 
power-law index a. Note that the vertical axis is represented 
by Ti5o- 1 like as Figure|8] As a reference, the 1D-HD model 
with F\ n /F K f = 1 is focused, where F K f is a reference value of 
the inflow energy flux which is equivalent to the model with 
Vjn/c= 1/10. This model provides the almost same physi- 
cal properties of the shock as those observed in the 2D-MHD 
model with a fixed shearing velocity 0.1c. This is verified by 
Figure [14] in which the result of the 2D-MHD model is pre- 
sented by the filled circle. 

The Ti5o-a relation in the 1D-HD model can reproduce that 
in the 2D-MHD model. The condition for the plasma velocity 
being relativistic in the 1D-HD model is also same as that in 
the 2D -MHD model, that is a > 5.5. These indicate that the 
shock driven outflow is accelerated to relativistic velocity by 
a purely hydrodynamic process even in the 2D-MHD model. 

Figure [15] demonstrates the r^o-a relation for the models 
with different inflow energies. The cross, filled triangle, and 
filled diamond denote the models with Fi n /F K { =1,1 /27, and 
1/1000 (that is, v in /c = 1/10, 1/30 and 1/100) respectively. 
The Lorentz factor becomes higher for the model with the 
larger inflow flux. It is important to stress that the condition 
for the plasma velocity being relativistic (Tiso > 2) changes 
depending on the amount of the inflow flux. When setting 
the lower value of the a for providing the relativistic velocity 
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FIG. 16. — The time evolution of the Lorentz factor in ID hydrodynamic 
simulation where a = 7.25 and v In = 0.1c. A solid line represents the time tra- 
jectory of the Lorentz factor at the shock surface measured in the laboratory 
frame. The meanings of vertical and horizontal axes are same as Fig. 9. 



as atrans, we can find a tra ns = 5.5, 6.0 and 7.0 for the model 
with Fi n /fref = 1, 1/27 and 1/1000 respectively. The transi- 
tion value atrans, at which the shock velocity undergoes tran- 
sition to the relativistic velocity, depends on the inflow flux. 
The larger inflow flux provides the smaller transition value. 

Note that the transition value atrans depends on the loca- 
tion where we measure the plasma velocity. We expect that 
the shock driven outflow is accelerated and finally reaches the 
relativistic velocity for all the continuous injection models as 
long as the initial density profile has a index a at least larger 
than 3 (the smallest value for our parameter survey). With 
fixed inflow flux, the transition value a^ans thus is expected 
to become smaller if the plasma velocity of the shock is mea- 
sured at a more distant location. In other words, transition 
value is a function of the inflow flux and the location where 
the plasma velocity is measured. 

Separately from the transition value, Figure Q3] indicates 
that there exists a critical value of the index a, at which the a- 
dependence of the plasma velocity changes. The critical value 
a cr it slightly shifts to the larger value with decreasing the in- 
flow flux, that is a cr it = 4.0, 4.25 and 4.75 for the models with 
Fi n /F re f =1, 1/27 and 1/1000 respectively. It is important to 
stress that the transition value a tra ns does not necessarily coin- 
cide with the critical value a cr i t according to the 1D-HD mod- 
els, even though it has a coincidental match with the critical 
value in the 2D-MHD models, (a^ans — a^rit — 5 see § 3.3.). 
The degeneracy of these two values would be dissolved when 
the input energy is varied in the 2D MHD models. 

To specify how the energy inflow pattern affects the dynam- 
ics, we examine the model with the impulsive inflow flux. The 
injected energy flux is fixed as Fi n /F re f = 1/1000 here. The 
open squares in Figure [TBlrepresent the 1D-HD models with 
the impulsive inflow. Note that the termination time, at which 
the energy injection is terminated, is fixed as lOr and is de- 
termined to provide the same total input energy as that in the 
2D-MHD model with the impulsive injection pattern (case B). 
We find that the different inflow pattern provides a different 
a-dependence of the Lorentz factor in the range a < a cr it- In 
contrast, there is little difference in the a-dependence of the 
Lorentz factor between two models when a > a cr if These 
indicate that the energy inflow continuously injected into the 
system does not affect the acceleration of the shock wave as 
long as the initial density profile is sufficiently steep. 



The evolution property of the shock is determined by the 
total amount of energy injected into the system in the dura- 
tion till the termination time in the impulsive inflow cases. 
By contrast, in the continuous inflow cases, the shock is con- 
tinuously accelerated without being terminated. The steady 
inflow cases thus provide the larger Lorentz factor than that 
for the impulsive inflow cases in the range a < a cr j t . 

Figure [16] shows the time evolution of the shock velocity 
when Fi n /F re f = 1 and a = 7.25 for the 1D-HD model. The 
axes have the same meanings as those of Figure [9] The dif- 
ferent curves denote the shock profiles at the different phases 
t /t = 80, 120, 160 and 200 respectively. The solid line repre- 
sents the trajectory of the Lorentz factor at the shock surface. 
The self-similar property of the relativistic shock in the 2D- 
MHD model is reproduced by the 1D-HD model. The self- 
similar evolution of the shock observed in the 1D-HD mod- 
els also follows the relation r sn oc r sn as is obtained from the 
2D-MHD models. These verify that a purely hydrodynamic 
process is responsible for the self-similar acceleration of the 
relativistic shock we found in the 2D MHD model. 

5. DISCUSSION 
5.1. Comparison with Other Self-Similar Solutions 

Piran et al. (1993) provides a detailed investigation of the 
dynamics of a fireball which is polluted by radiatively opaque 
plasma and whose internal energy is significantly greater than 
its rest mass energy. Because of its large opacity and huge in- 
ternal energy, the extremely hot plasma expands adiabatically 
as a perfect fluid. The homogeneous fireball accelerates ini- 
tially to a ultra-relativistic velocity and then coasts freely after 
converting all the internal energy into the kinetic energy. 

The Lorentz factor of the expanding fireball Tr, follows the 
relation Tfb oc r^ at the early acceleration stage where r^ is its 
radius (Piran et al. 1993). This relation is seemingly the same 
as the relation describing the forward shock obtained from our 
numerical models. In our model, however, the shock propa- 
gates not at ultra-relativistic velocity rather at mildly relativis- 
tic velocity (r srl ~ 10). In addition, unlike the fireball whose 
internal energy is significantly greater than the rest mass en- 
ergy, the internal energy of the plasma behind the shock is less 
than the rest mass energy in our calculations, as is shown in 
Figure[l7] Here we demonstrate the radial profile of the inter- 
nal energy and Lorentz factor of the plasma measured at the 
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FIG. 17. — The internal energy normalized by the rest mass energy and 
Lorentz factor of the plasma measured when t = 20()r at the equatorial plane 
in the model where a = 7.25 and @q = 3.2. The left and right vertical axes 
show the internal energy and Lorentz factor respectively. 



equatorial plane when t = 200r for the model where a = 7.25 
and fio = 3.2. The left vertical axis denotes the internal en- 
ergy normalized by the rest mass energy and the left axis is 
the Lorentz factor. The physics which characterizes the self- 
similar acceleration of the relativistic shock we found are thus 
essentially different from that of the fireball model. 

A spherical blast wave which expands supersonically into 
the surrounding medium is described well by a set of self- 
similar solutions. The Sedov-Taylor solution is the most fa- 
mous of these, but only describes the case where the velocity 
is in a non-relativistic regime (Sedov 1959). When the re- 
leased energy is much larger than the total rest mass energy 
of the explosion products and ambient gas, the velocity of the 
shock becomes relativistic. This type of shock also evolves 
self-similarly according to the Blandford-McKee solution at 
the ultra-relativistic limit (Blandford & McKee 1976). 

When we consider the shock driven by impulsive energy 
injection in a gas with the density profile p oc r~ a , it is decel- 
erated in the case a < 3 and accelerated for a > 3 (Shapiro 
1979; Waxman & Shvarts 1993). The evolution of the shock 
follows the second type self-similar solution when the density 
of the surrounding gas declines steeply with radius (Waxman 
& Shvarts 1993; Best & Sari 2000). In the ultra-relativistic 
case, the Lorentz factor of the shock T s h provides the relation 
T sh oc r m l 2 with m = (3 - 2^3)0; -4(5- 3 v^), where t is the 
time after the shock initiation (Best & Sari 2000). The index 
—m/2 ~ 0.77, 1 .0 and 1 .23 when the power law index a is 5, 
6 and 7 respectively. 

The physical conditions required for the solutions presented 
above are similar to ours. However, these are not accountable 
for the outflow property we found in our study. This might be 
because our self-similar solution is more suited for describing 
the mildly relativistic shock which can not be treated by either 
non-relativistic or ultra-relativistic solutions. We analyze our 
self-similar solution in more detail by comparing the other 
solutions in a separate paper . 

5.2. Application of Two-Component Outflow Model 

As an application of our two-component outflow model to 
the realistic astrophysical system, we consider the magnetic 
explosion event in the magnetar system. Giant flares from 
soft gamma-ray repeaters are believed to be powered by the 



explosive release of the magnetic energy on the isolated mag- 
netar (Thompson & Duncan 1995). It should be stressed that 
the giant flare on the magnetar is followed with an afterglow 
caused by the associated mass ejection event (Cameron et al. 
2005; Gaensler et al.2005; Taylor et al. 2005). 

Granot et al. (2006) diagnoses the observational properties 
of radio emitting ejecta observed in association with the giant 
flare on 2004 December 27 from SGR1806-20 (c.f., Taylor 
et al. 2005). It is suggested from their analysis that there 
is a possibility of the ejected matter containing two-types of 
components of sub-relativistic and mildly relativistic velocity. 
The origin and physics of the mass eruption event has not been 
determined although they give many clues to the model of the 
magnetar giant flare itself. Statistical analysis of giant flaring 
activity is needed for a better understanding of the magnetar. 

Although the numerical setting adopted in our model does 
not correctly capture the realistic features of the plasma sur- 
rounding the magnetar system, our model creates an out- 
flow with a mildly relativistic shock driven component and 
the subsequent sub-relativistic magnetically driven compo- 
nent loaded with dense plasma as a natural outcome of the 
magnetic explosion on a compact star. 

On the basis of our numerical model, the observed giant 
flare, energized from the magnetic explosion on the magne- 
tar, also powers the two-component outflow as is suggested 
by Granot et al. (2006). We aim to simulate the system of 
the compact star in the numerical setting more suited for the 
realistic magnetar system as part of our future work. 



6. SUMMARY 

We investigate the nonlinear dynamics of the outflow pow- 
ered by a magnetic explosion on a compact star using ax- 
isymmetric special relativistic MHD simulations. As an ini- 
tial setting, we set a compact star embedded in the hydro- 
static gaseous plasma which has a density p(r) oc r~ a and is 
threaded by a dipole magnetic field. A longitudinal shearing 
motion is assumed around the equatorial surface, generating 
the azimuthal component of the magnetic field. Subsequent 
evolution was simulated numerically. Our main findings are 
summarized as follows. 

1. Magnetic explosion on the surface of a compact star 
powers two-component outflow expanding into the static 
medium, i.e, magnetically driven and shock driven outflow. 
At the early evolutionary stage, the magnetic pressure due to 
the azimuthal component of the magnetic field initiates and 
drives a highly magnetized dense plasma outflow. Then the 
magnetically driven outflow excites a strong forward shock as 
the second outflow component. 

2. The expanding velocity of the magnetically driven out- 
flow depends on the Alfven velocity at equatorial plane on 
the surface of the compact star v A; o and follows a simple scal- 
ing relation v mag oc va,o 1,/2 - This scaling relation is accounted 
for in terms of coupling two balanced equations which are 
reduced from the induction equation and the energy equation. 

3. When the initial density profile declines steeply with 
radius, the outflow-driven shock is accelerated self-similarly 
to relativistic velocity ahead of the magnetically driven out- 
flow. The self-similar relation is T s h oc r s h, where T S h is the 
Lorentz factor of the fluid measured at the shock surface r s h. 
The pure hydrodynamic process is responsible for the accel- 
eration of the shock driven outflow and for the self-similarity 
of the shocked wave. 
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APPENDIX 

DERIVATION OF THE ENERGY CONSERVATION LAW ALONG THE STREAMLINE 

When the non-relativistic outflow is mainly polluted by the toroidal magnetic field, the radial outflow velocity should be 
controlled by the following equation of motion, in the spherical coordinate, 



/ dv r dv r 



-4(p + B -l 



B <p 2 G pM„ 
4-nr r 2 



(Al) 



where all symbols have their usual meanings and the final term of this equation represents a gravitational force of the central 
object with mass M*. Assuming the steady state, the radial velocity of the magnetically driven outflow should satisfy the steady 
state relation, which is reduced from equation (Al), 

\ i Tfl /r, 2 \ r , 2 ~\ 

(A2) 

By integrating this equation along the radial streamline from the equatorial surface /?„ to the arbitrary radius r and dropping the 
effect of the magnetic tension force on the dynamics, we can obtain the radial velocity v r (r) at the arbitrary radius r as follows; 
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When deriving this equation, it is assumed that the radial velocity at the stellar surface v r (r = R*) is negligible small (< 0.01c) in 
comparison with the outflow velocity v r (r) on the basis of our numerical results (see Figure 1), that is v r (R*) <C v r {r). 
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